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Abstract 



We present analytical and numerical study of discrete breathers identified as localized deformations 
of valence angles accompanied by change of valence bonds in crystalline polyethylene (PE). It is shown 
that such breathers can exist inside the optic wave number band and can propagate along the macro- 
molecule chain with subsonic velocities. Analytical results are confirmed by numerical simulation using 
Molecular Dynamics procedure. We examine also stability of the breathers relative to thermal excitations 
and mutual collisions as well as to collisions with acoustic non-topological and topological solitons. The 
conditions of breathers existence depend strongly on their generating frequency, relationship between 
stiffnesses of valence bonds and valence angles as well as on type of nonlinear characteristics of intra- 
chain interactions. 



1 Introduction 

In spite of growing interest to study of localized nonlinear oscillatory excitations, almost all studies in this 
field are devoted to straight chains. The reason is that even the most simple realistic models of strongly 
anisotropic systems, e.g. polymer crystals, turn out to be substantially more complicated systems than 
commonly considered one-dimensional models. 

Two-dimensional linear dynamics of planar zigzag chain was considered more than 60 years ago by Kirk- 
wood [1] in application to polyethylene macromolecule. From the other side zigzag and helix chains can be 
also considered as the discrete models of corrugated and spiral mechanical systems [2]. Growing interest to 
nonlinear dynamics of polymer chains and crystals in last decades is stimulated by numerous applications 
to such physical problems as dielectric relaxation [3, 4, 5], premelting and melting [6], heat conductivity 
[7, 8], polarization [9], fracture [10, 11], plastic deformation [12, 13, 14, 15], reading the information in 
DNA [16]. All these applications are based on two types of solitons : supersonic solitons of tension and 
compression [10, 17, 18] and combined soliton of tension (compression) and torsion [10, 19, 20, 21, 22], 
both types being in long wavelength region. 

Besides there are two examples of breathers in the models of polymer chains. 

First one is localized nonlinear vibration of C-H valence bond in carbon-hydrogen chains due to anhar- 
monicity of corresponding potential [23]. Such a breather has rather high frequency (~ 3100 cm" 1 ) because 
it is higher than lower optic branch of IR spectrum. Second example is localized periodic change of valence 
angles in PE chains accompanied by change of valence bonds [24]. This breather revealed by numerical 
simulation has significantly lower fundamental frequency (~ 800 cm" 1 ) that is close to lower boundary 
of first optic zone for PE crystals. However, one type of motion (namely transversal one) turned out to be 
predominant in the localization region of this breather because its spatial frequency is close to boundary 
value of wave number (k = ix). Besides, only stationary breathers have been revealed. Meanwhile, the 
valences angle potential for PE chain presented in [20] demonstrates the case when lower boundary of optic 
dispersion curve is inside of wave number diapason and is not close to its boundaries. It means that possible 
breather in this case has not a predominant component (longitudinal or transversal) and its study becomes 
much more complicated. 

In this paper we remove both restrictions mentioned above. Besides, we present first analytical descrip- 
tion of short wavelength nonlinear dynamics of PE chain in all diapason of optic frequencies alongside with 
numerical investigation of breathers stability relative to thermal excitations as well as mutual collisions and 
collisions ith acoustic non-topological and topological solitons. It worth mentioning that optic breathers 
in PE crystal weakly depend on interchain interaction (contrary to topological solitons). However, this 
interaction was taken into account in the course of numerical studies. 

The commonly studied breathers in attenuation band can be, in principle, linearly stable [25], [27]. How- 
ever as it was shown in [24], their nonlinear interaction with extended modes (phonons) leads to finite times 
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of life in realistic models of polymer crystals for breathers. Therefore, in such models there is not quali- 
tative difference between the breathers in attenuation and propagation bands, if naturally they can exist in 
propagation band at all. From other side, the finite time of life does not mean that such breathers are not 
of physical importance. They can be manifested in different physical properties [26] and their instability 
relative to interaction with phonons leads to non-monotonous dependence of breathers contribution with 
increasing the temperature [24]. 



2 Description of the model 

Because in numerical study we deal not only with planar motion of the chain but take also into account the 

spatial movement and interchain interactions, three-dimensional hamiltonian is presented below. 

It is suggested that the zig-zag initial structure of PE macromolecule is directed along the ,2-axis in plane 

(x, z), while y-axis is perpendicular to the plane (x, z) so that (x, y, z) is the positive global reference of the 

problem, 8 Q being the angle of the zigzag chain. 

Initial coordinates of the n-th mass of the chain are : 

z° = nl z , x° n = (-l) n+l l x /2, (1) 

where l z = sp and l x = cp [s = sin(6> /2), c = cos(6> /2)) are longitudinal and transversal dimensions 
of the initial zigzag chain. We use approximation of "united atoms" because relative motions of hydrogen 
atoms are non- significant when dealing with backbone deformation predominantly. So the magnitude of 
every mass is supposed to be equal to 14 a.u. It is convenient to introduce the relative coordinates 

%n ^ni % n Z n , (2) 

where w n , u n , v n determine the longitudinal, transversal and out-of-plane displacements of the n-th mass 
from its equilibrium state. Let us denote by p n and 9 n the current length of the valence bond and current 
valence angle respectively. We introduce also <p n as the angle between d n _i and the plane generated by d n 
and d n+ i (conformation angle). 

The dynamics of the chain is governed by the following hamiltonian function of the chain [24]: 

+00 ^ 

n=—oo 

+ [Pl + P2 COS(0„) + p 3 COS(30„,)] 

+-A[cos(0 n ) - cos(# )] 2 

+V[1 - e-i^-ri] 2 } + Z(u n , v n , w n ) (3) 

where dots denote the time derivation, A and V correspond to energies of valence angles and valence 
bonds respectively, and the last term the energy of interaction of the n th unit with six neighboring chains 
(substrate potential). Parameters Pi,P2,P3 satisfy conditions p 1 = p 2 + £>3, P21P3 > 0, so that undeformed 
state corresponds to 4> n = n. Let us set ki = 2Vq 2 and k 2 = A sin 2 9 . Then 

Pn = ( U n+1 - Un- + On+1 ~ V nf + On+1 ~ W n + l z f ■ (4) 

We have also: 

COS# n = (— a n -l,l«n,l — an-l,2fln,2 

—0'n-l,2,0"n,'i) / PnPn~l, (5) 

cos0„ = {b nA b n+1A + b nt2 b n+h2 

+b nj3 b n+lt3 )/fl n fl n+1 , (6) 
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where 

b n ,l — a n-l,2 a n,3 ~ Qn,2°n-1,3) 

b n ,2 = On-l,30 n ,l — «n,3 a n-l,l, 

b n ,3 = 0<n-l,X a n,2 — «n,l«n-l,2 

are related to inner product vector coordinates, and (3 n = (b^ 1 + b 2 n2 + 6 2 3 ) 1//2 are norms of corresponding 
vectors with 



a n,l 


= Un+1 




O n ,2 


= V n+ i 


- v n , 


a n,3 


= w n+1 


-w n + l z , 


Pn 


= («n,l 


+ <2 + < 3 ) 1/2 



which are the vector coordinates. The substrate potential 

7TU) 71 IV 71 IV 

Z(u, v, w) = e w sm 2 (nw/l z ) + -K u [l+e u sin 2 (—)]{u--l x [l-cos((—)]} 2 + -K v [l + e v sin 2 (— )]t; 2 , 

where e u = 0.0674265 kJ/mol, e v = 0.0418353 kJ/mol, e w = 0.1490124 kJ/mol, K u = 2.169513 kJ/ A 
mol 2 , K v = 13.683865 kJ/ A mol 2 . 

3 Planar motion of zigzag chain 

We consider further analytically only in-plane dynamics of zigzag. In both linear and nonlinear problems 
such a dynamics can be fully separated from out-of-plane motion. Therefore, general Hamiltonian can be 
simplified and one can obtain in both physically and geometrically nonlinear approach: 

+oo 

H = ^2 -m(u 2 n + u> 2 ) terms up to order 4 from 

n=— oo 
f + °° 1 

£ -A[cos(0 n ) - cos(^ )] 2 

ln=— oo 

+oo ~| 

+ ]T V[l-e- q{pn - po) f\ (7) 

n=— oo J 

There are two systems of parameters describing PE macromolecules which are supposed in [19] and [20] 
respectively which differ by value of parameter k 2 (k 2 / sin 2 9 = 130.1 kJ/mol in [19] and 529 kJ/mol 
in [19]). We will consider both systems supposing [18, 19] that ki/q 2 = 334.7 kJ/mol, q = 19.1 nm™ 1 , 
6 = 113°, po = 1-54 A so that 5 = k 2 /k lP 2 = 0.019 [19] or 0.078 [20]. 

The deformations are presented by their power expansions including the terms of first and second order 
with respect to displacements: 

Ap n = Pn~ P0 = 

(u n - u n+l )c + (w n+1 - w n )s + — ((«„ - w„+i)c + (w n+1 - w n )sf + . . . , (8) 

^Po 

= n — 6q = 

s c 
— (2u n - n n „i - u n+l ) H O n _i - w n+1 ) 

Po Po 

H 2 [{u n ~~ u n+lf + (U n - U n _i) 2 - (w n - W n+ i) 2 - (W n - W n ^i) 2 ] 

Po 

C 2 - S 2 

H 2 \-( Wn ~ W n -i)(u n - U n -!) + (W n+1 - W n )(u n - U n+1 )} + ... (9) 

Po 
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Equations of motion are obtain in the form: 



d 2 dH d 2 dH 

Ot OUr, Ot l OWr, 



for the nth particle. 

Linearized equations are given by 



d 2 

m—u n + kxsciwn+i - w n _i) + kic 2 {2u n - w n _i - u n+1 ) 



8t 2 



+ k 2 cs(w n _ 2 + 2w„_i - 2w n+1 - w n+2 )/p 2 

+ k 2 s 2 (u n - 2 - 4Mn_l + 6u n - 4u n+1 + u n+2 )/pl = (11) 



d 2 

m ~^ w n - hs 2 (w n+1 - 2w n + + kisc{u n+ i - U n _i) 



dt 2 



- k 2 c 2 (w n - 2 - 2w n + w n+2 )/pl 

+ k 2 s 2 (-u n „ 2 + 2u n _i - 2u n+1 + u n+2 )/pl = (12) 



4 Dispersion relations 

The linear dispersion curves have well known form (see, e.g. [24]), which is described by relation : 



U 2 (k) = f{k) = UJ 2 (k) ± yju&k) ~ Uf(k) (13) 

where 

uJo(k) = Ci(l — cos^o cosk) 

+2C 2 (l-cosk)(l + coskcos6 ), (14) 
ujf(k) = 8CiC 2 (l - cos A;) sin 2 A;. 

Here, Ci = ki/m, C 2 = k 2 jmp\, k = kl z , k is wave number. Signs "minus" and "plus" correspond 
to acoustic and optic branches respectively. We consider further the optic branch only. Corresponding 
asymptotic representations of the linear frequencies can be presented as follows in the vicinity of arbitrary 
wave number k*\ 

u 2 (k) = n 2 + u(k - it*) + p(k - k*) 2 + ... = f(k), (15) 



with 



^9 ^ r ,„ s ,,« n 4C 2 cos 2 (k*) sin 2 (9 ) . 

» = 2C '" - "»W "»(*')! ' i 1 cos(W cos(l) " - «"(* "■ 
i/ = 2Cx sin(fc*) cos(^o) 

+ 2C 2 sin(2A:*) sin 2 (6>o) [(1 _ cQs(r))(2 _ ^ + _ ^ 

(1 — cos(#o) cos^/c*))^ 

/x = d cos(flo) cos(fc*) + Sm2( y [-4 + 12 cos(fc*) 

( cos (t^o ) cos(fc*) — l) 6 

+4cos 2 (A;*)(2-3cos(^o)) -2cos 3 (A;*)(-2cos 2 (^o) +3cos(^ ) +9) 
+2cos 4 (fc*) cos(fl )(cos(#o) + 11) - 8cos 5 (^)cos 2 (^ )]- 
Table 1 provides values of Q 2 , v, p for some particular k. 

Analysis of general expressions of dispersion curves in application to realistic values of chain parameters 
reveals two types of behavior. The corresponding plots of dispersion curves are presented in Fig. 2. 
Let us introduce C 2 = 5C\. The signs of numerical values of optic dispersion curves curvature H are: 
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Table 1: Values of f2 2 , u, fi for some particular k. 



k 


n 2 


V 


/' 










Ci(cos(# ) +4c 2 5) 


TV 


4C 1 (c 2 + 4s 2 <5) 





cos(^ )c 2 + 20cV<S 
1 c 2 + 4s 2 5 


tt/2 


2d 


2d cos(fl ) 


4Cl %l-2C 2 



• 5 = 0.019, H(Q) < 0, W(tt) > 0, H(ir/2) < 0, 

• 5 = 0.078, W(0) < 0, H(tt) < 0, W(7r/2) > 0. 

We can see that for k ~ curvature is negative in both cases. For A;^7r, A;^7r/2it could be both positive 
or negative depending on the magnitude of 5. The signs of the curvature determine as we will see the signs 
of second spatial derivatives in final continuum equations of motion. Let us define k G]0,7r[ so that the 
corresponding coefficient v is equal to zero: f'(k ) = (k is defined only for 5 = 0.078, because for 
5 = 0.019 the value v = at boundaries of wave number diapason only). 

5 Introduction of modulating functions 

In this section we introduce slow modulating functions in order to extend analysis based on normal modes 
that can be done for linearized equations to the nonlinear case. So we consider equations for u n , w n both 
with modulating continuous space functions W, W expanded close to arbitrary n th particle: 

/ ,^rr 9U 1 2 2^1 / , N rfr 1 2 2 8 2 U . 

u n+m = cos{mk)[U + me -Q£ + 2 m e "^2"J + sm(mk)[U + + -m e -^\ + ... (17) 

, ,, rTJ7 dW 1 2 2 d 2 W, , ,, r ~ dW 1 2 2 d 2 W, 
Wn+m = cos{mk)[W + me_ ^- + -m e -—-J + sm(mk)[W + me ~^~ + -m e -^-] + . . . (18) 

where e is a small parameter characterizing distance between particles in the units £ = s~ l z/l z , z being 
dimensional distance. 

In the same spirit, we consider equations for n n+ i> Wn+i> both with the modulating functions expanded 
close to the n + 1 th particle (i.e. now with u n+ i = U, w n+ i = W): 

< ^nr dU 1 2 2< 92f/ l , rMfr d ^ 1 2 2^1 

u n+1+m = sm(-mk)[U + me ^ + -m e -^-\ + cos(mk)[U + m ^-g^ + -m e -^-\ + ... (19) 

INrTTr dW 1 22 d 2 W, , dW 1 2 2 d 2 W, 

w n+1+m = sm(-mk)[W + me_ ^- + -m e -—-J + cos(mfc)[H/ + me_ ^- + -m e -^-] + . . . (20) 
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Necessity to consider the expansions in the vicinities of two points is caused by twice degeneration of linear 
normal normal modes spectrum (except the boundary values of wave number). By substituting modulating 
functions in the linearized equations for u n , u n+1 , w n , w n+1 (11), (12), one obtains 4 partial differential 
equations in the general case that take into account linear part of the equations for modulations of normal 
modes. Linearized equations for modulations of normal modes have to be in accordance with dispersion 
relations presented above. We present them below for special case k = ir/2. 



d 2 W 

+sc(2C 2 + C^e 2 —— + 2(c 2 d + 2s 2 C 2 )U + 2sc(2C 2 - C X )W = 0, 
9 2 w n „ 2 dw A „ dU A „ 2 2 d 2 W 

i^- 2C ^ +4csC2e ^ +4C2Ce w 

+sc(2C 2 - C x )e 2 ^r + 2(s 2 d + 2c 2 C 2 )W + 2sc(2C 2 - d)U = 0, 
d 2 W 

+sc(2C 2 - C^e 2 —— + 2(c 2 d + 2s 2 C 2 )U + 2sc(2C 2 - C X )W = 0, 

^ + 2C lS 2 e—-AcsC 2 e- + AC 2 c 2 e 2 w 
+sc(2C 2 - C^e 2 ^ + 2(s 2 C 1 + 2c 2 C 2 )W + 2sc(2C 2 - C X )U = 0. 



(21) 



(22) 



(23) 



(24) 



d 2 d 2 

Then we introduce operator approximation for -^—;U n and -—ru n+ i in corresponding equations using previ- 

ot at 1 

ous Taylor expansions and expansion of dispersion relation according to 

9 2 2rr r,2rr ^ 2^U 

-u n ^-. 2 U = -n 2 U + ue- + ,e 2 w ... 

92 2fr r»2rr 9U K > 

l -u n+1 ^-uU = -nU-ve— + f ,e— + ... 

Substitution of (23), (24) into starting nonlinear equation (9) leads to four nonlinear partial differential 
equations (NPDE) with respect to functions W, W, U, U. To reduce the order of this system one can 
express the functions U, U via W, W using the equations for U and U, and taking into account the terms 
up to third degree. As this takes place, we have to replace the nonlinear terms with respect U, U by their 
expressions via W, W obtained from (21), (23), (24). Then NPDE for U and U provide two independent 
linear relations relative these functions that can be solved to obtain U and U versus W and W. Indeed 
truncation of the full NPDE at order 3 are used. At this step, one can notice that for k = ix /2 or k = k 
displacements u n and w n are of same order: we have longitudinal and transversal motions. If measuring 
the length in units l z we have k ~ %/2 : u n ~ w n . In strongly coupled general case, the relations can be 
written as: 

dW dW od 2 W nd 2 W 
U = A x W + A 2 W+F l2 {W, W)+F 13 (W, W)+A 3 6— + A 4 6— + A 5 6 2 — + A 6 e 2 — + .. . , (26) 

dW dW d 2 W ^d 2 W 
U^B 1 W+B 2 W+T 22 (W,W)+T 23 (W 1 W)+B 3 e— + B i e— + B 5 e 2 — + B 6 e 2 — + ..., (27) 
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with Tij, i = 1, 2 of degree j in W, W. Ak, Bk, k = 1, . . . , 6 are given constants. For example, if k — ir/2, 
one has: 

f Ai = 0,A 3 = — ,F 12 (W,W) = D(W 2 + W 2 ), 
s 



B 2 = 0,B 1 



D 



,T 22 (W,W) = -D(W 2 + W' 
3cqp C! 



(28) 



2 S 4(d - 2C 2 ) 

For = or k = ir, degeneracy of order occurs: u or w displacements dominate. We have: 

k ~ : u n ~ fciu n ~ 5w n , we suppose that k ~ e « 1 

/c ~ 71 : U! n ~ (7T — /c)« n ~ £W n . 

Then modulating functions introduced in equations for u n and w n provide only 2 NPDE. In these two 
particular cases, one needs only to consider equation for u n (or respectively w n ). The relations, connecting 
u n and Wn are obtained in the form 

~cdW 



u n = U = e 



s d£ p s 



+ -(2 - 3c 2 + 3s 2 p q - 85s 2 )W 



+ ... 



for k = and 



for k = i\. 



w n = W 



sc 1-45 dU 



2 A5s 2 + c 2 d£ 



(29) 
(30) 



6 Nonlinear partial equations 

Let us introduce new non-dimensional time by relation r = Qt, where f2 is linear frequency for consid- 
ered normal mode. After substitution of the expressions of the displacements via modulating functions 
into equations for w n and w n+ i, one can obtain corresponding nonlinear partial differential equations for 
modulating functions in the form: 



d 2 W f'(k) dW f"(k) 2 d 2 W 



cV 2 



n 2 d£ n 2 a? 



(3D 



+—[D 1 W 2 + D 2 WW + D 3 W 2 + D 4 W 3 + D 5 W 2 W + D 6 WW 2 + D 7 W S 



tt 2 



d 2 W f(k) dW f"(k) 2 d 2 W ~ 



1 

"Q 2 



dr 2 tt 2 d£ 
D 2 WW - D X W A + D 7 W 



—D?.W 2 n T;r/T;r/ n T;r/2 1 n T;r/3 



a 2 a? 

D 6 W 2 W + D 5 WW 2 



(32) 



D 4 W £ 



0. 



where Dj,j = 1, ... ,7 are constants depending on C±, C 2 , po, q, 6q. For k = n/2, the expressions are given 
in table 2. They are not given for the general case because of size of expressions. Numerical values for 
K = K , in re-scaled nondimensional form corresponding to e = 0.05 are given in table 2. 

Again, degeneracy of cases k = or k = ix leads to only one equation obtained from equation for w n 

as: 

d 2 W d 2 W 

Ai£ 2 ^-r + W+a 2 W 3 



dr 2 



Ai 



C x cos(^o) + 4c 2 C 2 



0. 



4s 2 d 



(33) 



a 2 



1 ro/c 2 /)\2/2 i 14 222 

[2(5c - 4)c /s + yp g s 



12p gc 2 + 325c 2 
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Table 2: Values of constants Dj, j = 1, . . . , 7 for k — 7r/2 and numerical values for k = k , corresponding 
toe = 0.05. 



j 


k = tt/2 


k = k 


1 


(d - 2C 2 )Dcs/C 1 


2.25 


2 





-5.74 


3 


(Ci - 2C 2 )Dcs/C l 


2.57 


4 


-D(3C 1 qp + C 1 -2C 2 )/C 1 


1.89 


5 


D(3C iqPo + d- 2C 2 )/C 1 


-1.14 


6 


-D{3C l q P o + C l -2C 2 )/C l 


2.42 


7 


D{3C iqPo + Ci- 2C 2 )jC x 


-1.14 



for k = and 



A 2 e 2 ^r + U + a,U 2 + a 2 U 3 



dr 2 



d cos{6 )c 2 + 20c 2 s 2 C 2 
'' ~ 4(c 2 + As 2 5) 2 ' 

3c 2 (s 2 -c 2 qp )+As 2 (l-Ac 2 )5 
c c 2 + 4s 2 S 



(34) 



a 2 



2[3cV(5c 2 - 1) + 18cVgp - 7c 6 g 2 pg - 12£s 2 (48c 4 - 24c 2 + 1)] 

3c 2 (c 2 + 4s 2 <5) 



for k = ix. In the last NPDE, nonlinear terms involving spatial derivatives can be neglected because they 
are of higher orders by e. 

All these NPDE have been given in dimensional form. They could be re-scaled according to physical 
orders of the different variables (e.g. replacing the variables U, U, W, W by epoU, epoU, epoW , epoW). 
This has been done for numerical applications in the cases k « 0, k « ix and k = ko below. Such re-scaling 
makes the orders of nonlinear terms to be similar to those of the terms containing second spatial derivatives. 



7 Transition to complex variables 

The NPDE of previous section can be written under first order by time if setting 

fc~0, V(£,t) = (W' + iW), **(£,t) = (W'-iW), 

k~n, *{Z, T ) = {U' + iU), **(£,T) = (U'-iU), 

*i(£, t) = (W> + iW), r) = (W - iW), 



k ~ tt/2 or k , 



We obtain two NPDE, for k = 0: 



tt 2 (£, t) = (W' + iW), t) = (W' - iW). 



Ax 2 d 2 



Conjugate equation , 



(35) 



(36) 



9 



and k = 7T 



= - -e 2 |^(* - **) - - #*) 2 _ ^ _ ^)3 + 



r .0^ 
97 



2 

Conjugate equation , 
In general case one can obtain four NDPE, e.g. for A; = n/2 they have linear part 

r Mi T /'(A;) a 2 T ., f n (k) 2 d 2 /T 

9r 1+ ft 2 2j ^ 2 lj ' 

Conjugate equation, 

M 2 T /'(*0 <9 2 , T T ., /"(A;) 2 9 2 /T 
i = - - e f^i - \EC) - — — e f\Po - 

Conjugate equation. 
Let us introduce the changes of variables 

T =T, T\ — ET , T 2 = E 2 T , . . . 



(37) 



(38) 



(39) 



Setting ij) = cpe lT (for k = 0, 7r) and = 0-,-e"", j — 1,2 (for k = %/2) we use farther the power expansions 



e( (po + e<pi + £ 02 



and respectively 



0j = £ ( ( l ) jO + + ^ 2 0j2 + • • •) 

Then we obtain at order e° the relations: 



(40) 
(41) 



k = : -% 



k = ir : — z 



<9r n 



<9r 



0. 



(42) 

(43) 
(44) 

<9r <9t 

It means that for k = 0, ir, O = O (£, r 1; r 2 ) and for k = k , 1)O = 0i, o (£, n, r 2 ), 2 , o = 02,o(£, n, r 2 ). 
For k = 7r/2 presence of first space derivatives lead to different problem. 
At order e l we have the equations 



k = ko : — i - 10 =0, — i - 20 = 0, 









"<P0 






901,1 


<90i,o 


9r 


9n 



0. 



'2,1 . U(p 2 fl 



dr dr t 

Conditions of absence of resonances lead to equations 



k = 0, 7i : 

dry 







0. 



(45) 

(46) 
(47) 

(48) 
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k = k : ^ = 0,^ = 0. (49) 

OTi OTi 

Therefore O = 0o(£, r 2 , r 3 , . . .) for k = 0,n as for the case k = k , we can write that 0j >o = 
0i,o(£,T 2 , . . .). j = 1,2. 

At order e 2 the conditions of absence of resonance terms lead to final complex equations in main ap- 
proach: 

/c = 0,7r: -«^r + P^- + « I 0o I <Po = 0, 

k = 0:a = -^,P = \ 1 /2, (50) 
8 

3a 2 5a? „ , . 
fc = 7r,a = -^ + ^,/3 = A 2 /2. 

which is Nonlinear Schrodinger Equation (NSE) for k = and k = it respectively. For k = k , the resonant 
equations have the general form: 

■ <901,Q A 2 (9 2 0iq 2 2 2 

— + -Pi | 010 | 010 + P2| 010 I 020 + P5| 020 | 010 ^ 

+Pel 020 | 2 02O + P40?n02O + P509O01O = 0, 



A 2 d 2 - 



. 0^2,0 / >2 U (P20 n I , ,2, D , , ,2, D , , ,2, 

~ Z ~^7 ^ + Pi | 020 | 020 + P2| 020 | 010 + P5 1 010 | 020 ^ 

+P 6 | 010 iVlO + P 4 0l O 0^ + P50 2 O02^ =0 = 0. 

with Pi, i = 1,6 constants derived from NPDE (31) and (32): 

Pi = -D 4 --D 2 1 + -D 2 D 3 , 
8 12 1 24 3 ' 

P 2 = -D 5 + -D 2 -—D 1 D 2 + —D 2 
8 2 3 12 12 2 

P 3 = — D 2 -— D 1 D 2 + -D% + -D 5 , 
3 12 3 24 8 2 8 

^ = -IdI + ^d.d, + ~D 2 D 3 + ±D 6 , 

p 5 = Id 7 -^d 2 d 3 + ^ Di d 3 , 

P 6 = -^i^s - ^f>\ + \d 6 + l -D 2 D 3 + -^D X D 2 . 



8 Soliton-like solutions 

As it is known, the type of soliton-like (breather) solutions depends strongly on relationship between the 
signs of constants a and (3 of NSE. When these signs are similar, NSE admits the envelope solitons. In 
opposite case, NSE admits "dark" solitons. From this point of view there is a significant difference between 
two systems of zigzag parameters introduced above. 

For the case k ~ 0, and both values of 5, the signs of a and (3 are the same (< 0) (see numerical values 
in table 3). 

For the case (k ~ ir, 5 = 0.019) the dispersion curve has a form similar to Figure 2 (a). In such a case 
the signs of a and (3 are the same (> 0). Then one can obtain the the envelope solitons (or breathers) as 
particular solutions: 

o (£,r 2 ) = (2A/a s ) 1/2 exp(z</2 v /^ + 

xsech[A 1 / 2 (e/ v /^ + ^r 2 )], (53) 
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Table 3: Values of a, (3 versus 5 for k ~ or it. 



k 


5 


a 


f3 


o 


0.019 


0.053 


-0.07 





0.078 


0.054 


-0.05 


71 


0.019 


0.035 


0.021 


71 


0.078 


0.033 


-0.23 



Table 4: Values of a, (3 versus 5 for k ~ k , 7 = ±1. 



7 


5 


a 


P 


+ 1 


0.078 


-0.010 


0.28 


-1 


0.078 


2.47 


0.28 



where uj = v 2 /4 — A. 

For the second system (k ~ 7r, 5 = 0.078) the dispersion curve has a form similar to Figure 2 (b) (optic 
branch) and dark solitons exist. 

In the case k ~ ko, breather or dark soliton can exist near minimal frequency. This fact has been 
confirmed by computer simulation. Looking for particular solutions of the form 2 ,o = 70io> we obtain 
only two possible values for 7: 7 = ±1. The two coupled complex NPDE (51) and (52) can be reduced to 
only one Schrodinger equation of the form 

~ % -^- +y-Q£T + a \ 0io I 0io = 0, (54) 

with j3 = — A 2 /2, a = Px + P 4 + P 6 + ^f(P 2 + P3 + P5). Numerical values are given in table 4. The breathers 
exist for 7 = 1 and dark solitons for 7 = — 1. 



9 Numerical simulations 

To check the validity of assumption made in the analytical study, we have undertaken a numerical treatment 
of the breathers existence as well as their stability in free motions, under collisions and thermal perturba- 
tions. 

While numerical modeling the breathers and their dynamics, we consider the following system of equa- 
tions corresponding to Hamiltonian (3): 

d 2 dH d 2 dH d 2 dH 

ot ou n at 1 ov n at 1 ow n 

forra = 1,2,..., N. 

We use initial conditions in agreement with approximate analytical solution. Because of small difference 
from exact solution, there will be a phonon radiation. For its absorption, the viscous friction is introduced 
at the end of the chain. As it was mentioned above, we deal with two systems of parameters for PE crystal. 
One of them is[19]: 

m = 14m p , pi = 8.37kJ/mol, p 2 = 1.675kJ/mol, 
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p 3 = 6.695kJ/mol, A = 130.122kJ/mol, (56) 
00 = 113°, D = 334.72kJ/mol, 
q = 1.9lA~\ p = 1-53A, 

where m p is the proton mass. Second system of parameters which is used in [20] differs with more high 
value of parameter 

A = 529kJ/mol. (57) 

When using the first system of parameters (56), a geometric nonlinearity plays a crucial role in nonlinear 
dynamics of the PE chain. In the second case (57) a physical nonlinearity becomes more essential. More- 
over, the dispersion curves for these two cases have different view. Respectively, the optic breathers will 
have different view and different regions of existence in parametric space (Fig. 2). Therefore we study 
numerically their properties for both systems of parameters. 

Numerical integration of the equation of motion (55) has confirmed that in accordance with analytical 
study, the optic breathers exist for both systems of parameters (56) and (57) near low boundaries of the 
frequencies of optic phonons. Typical distribution of relative displacements in the localization region of 
planar breather is presented in Fig. 3, 4, 5. One can see that shift of frequency (see Fig. 3 and 4) leads to 
narrowing of breather. Analogous effect is achieved also when increasing the intensity of excitation - see 
Fig. 5. Characteristics of breathers for the systems of parameters (56), (57). are essentially different - see 
Fig. 6 and 7. In the localization regions of breathers the local change of valence angles is accompanied by 
local extension of zigzag (average values of relative longitudinal displacements u> n +i — w n are positive). 

To consider an interaction of breathers with thermal vibrations of the chain, N segments of the chain 
were inserted from every side into thermal bath with temperature T. Then the Langevin equations of motion 

d 2 dH . 

d 2 dH 

m—v n = h Vn -T n mu n , (58) 

dt 2 dv n 

d 2 dH > „ 

dt 2 dw n 

where the Hamiltonian of the system H is given by Eq. (3), £ n , r) n , and ( n are random normally distributed 
forces describing the interaction of nth molecule with a thermal bath, coefficient of friction T n = 0, forces 
£n, Vn, Cn = for A" < n < N - N and T n = T for n < N and N — N < n < N. Coefficient of friction 
r = l/t r , where t r is the relaxation of the velocity of the molecule. The random forces £ n , rj n , and ( n have 
correlation functions 

(UhMh)) = (vn(hMt 2 )) = (C„(*i)6(*2)) 

= 2mTk B T5 n i5(ti - 1 2 ), 

(UtiMh)) = (Uh)Ci{h)) = (r)n(ti)(i(t2)) = 0, 
1 < n, I < N , N - N < n,l < N, 

where ks is Boltzmann's constant and T is the temperature of heat bath. 

The system (58) was integrated numerically by the standard forth-order Runge-Kutta method with a 
constant step of integration At. Numerically, the 5-function was represented as 6(t) = for |t| > At/ 2 and 
5(t) = 1/ At for \t\ < At/2, i.e. the step of numerical integration corresponded to the correlation time of 
the random force. In order to use the Langevin equation, it is necessary that At t r . Therefore, we chose 
At = 0.001 ps and relaxation time t r = 0.1 ps. 

To avoid an effect of friction coefficient on the behavior of breather, it was isolated from heat bath. For 
this, the stationary breather was situated at the center of the chain with N = 50. In such a case, the breather 
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Table 5: Dependence of breathers time of life t a on the chain temperature T (frequency uj = 820.5 cm , 
parameter A = 130.122 kJ/mol and u = 928.4 cm" 1 , A = 529 kJ/mol). 



A (kJ/mol) 


T(K) 


1 


2 


3 


5 


10 


20 


130.122 


t a (PS) 


180 


94 


68 


46 


33 


22 


529 


t a (PS) 


1052 


625 


386 


228 


121 


68 



can interact only with thermal phonons arising at the ends of the chain, which are connected with heat bath. 
The numerical integration of the equations (58) has shown that, contrary to isolated chain, the breathers in 
thermalized chain have a finite time of life. However, this time is large enough to provide a significant role 
of the breathers in different physical processes. Breaking of the breathers in thermalized chain is shown at 
Fig. 8 and 15 (a). 

To estimate the breather time of life in thermalized chain let consider the temporal dependence of the 
energy (Fig. 9 and 13). It is seen, that the energy decreases monotonously by exponential low E{t) = 
E(0) exp(-at). One can determine the breather time of life as half of breaking period t a = ln2/ab 
(E(t a ) = E(0)/2). Dependence of t a on the chain temperature T is presented in table 5. We conclude, that 
breather time of life is proportional to ratio of its energy to temperature. 

Besides the considered breathers, the supersonic acoustic solitons can exist in PE chain [17]. Therefore, 
it is desirable to study their interaction with optic breathers [we choose further for definiteness the system 
of parameters (56)]. In such a case the acoustic solitons are caused by a local longitudinal compression of 
the chain (the properties of acoustic soliton were presented in [17]). 

The collision of acoustic soliton with the stationary breather is shown at Fig. 10 (a). It is seen, that 
considered interaction does not lead to noticeable change of soliton energy, but the breather acquires a 
small momentum in the direction of the soliton motion. E.g., after collision of the soliton with velocity 
s = 1.024so, where s = 7790 m/s - velocity of long wavelength optic phonons, the breather with frequency 
uj = 820.5 cm" 1 becomes to move with constant velocity s = 0.0255so- After every new collision, the 
breather velocity increases attaining the value s = 0.106so in limit. 

The motion of the breather is accompanied by phonon irradiation with low intensity that leads to de- 
creasing of breather energy. The temporal dependence of the breather energy is shown at 1 1. It is well seen 
that change of velocity and energy decreasing may be fixed at very large time only (~ 10 ns). 

It is necessary to note that propagating breather has the frequency slightly exceeding the low boundary 
of frequency spectrum. So, we deal here with the breathers in the propagation zone of optic spectrum. 
However, the time of life in this case turns out also to be large enough as well as in the case of breather in 
the gap between acoustic and optic branches of the dispersion curve. Actually, let us consider the collision 
of propagating breather (u = 830 cm" 1 ) with stationary one (u = 820 cm -1 ) [see Fig. 10 (b)] and collision 
of two propagating breathers {uj = 830 cm" 1 ) [see Fig. 10 (c)]. One can see that the breathers with 
frequency in propagation zone can move freely without any noticeable changes and interact with similar or 
stationary ones as elastic particles - they exchange by momentum without any energy decreasing. 

In thermalized chain both propagating and stationary breathers have the time of life, which is propor- 
tional to ratio of breathers energy to temperature and does not noticeably depend on the breather velocity 
(Fig. 12). 

We have shown that the optic breathers can exist in both isolated chain and chain interacting with 
neighbor ones in the crystals. As this takes place the characteristics of the breather are not depend noticeably 
on the intermolecular interactions. In this case we can use approximation of immovable neighbor chains. 
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Corresponding Hamiltonian (3) has a view 



+00 1 

h = E U m K + vl + <\ 

n=—oo 

±[Pi + P2 cos((p n ) + p 3 cos(3(p n )] (59) 

+ ^A[cos(# n ) - cos(^ )] 2 

+V[1 - e-ffC**-")] 2 + iy( Mn , Unj «,„)}, 

where the function W{u n , v n , w n ) describes the interchain interaction. Detail description of this approach 
is given in [21]. 

If taking into account interchain interaction, three types of topological solitons with topological charges 
q = (qi, q 2 ) appear. Their properties are described in [21]. The solutions of the first type have topological 
charge q= (±1,0) and describe a localized longitudinal extension (compression) of the chain by one period. 
The solutions of the second type have the charge q= (±0.5, 0.5) and describe a longitudinal extension 
(compression) of the chain by halved period and simulations twist by 180°. The solutions of the third type 
have the charge q= (0, 1) and describe the twist of the chain by 360°. All these topological solutions have 
subsonic spectra of velocities and can move without phonon irradiation. 

We considered an interaction of stationary optic breather having frequency u = 820.5 cm" 1 with mov- 
ing topological solitons. As it is seen from Fig. 13 such an interaction is elastic one if first component of 
topological charge q x < 0, i.e. if a compression in the localization region is absent. In opposite case, when 
qi < 0, the interaction with topological soliton leads to breaking the breather - see Fig. 14. Moreover, the 
bound state of breather and topological soliton (qi > 0) may be energetically profitable. Such a coupling 
increases the lifetime of breather in thermalized chain [compare Fig. 15 (a) and (b)]. 

10 Conclusion 

The optic breathers which are localized coupled longitudinal-transversal nonlinear excitations can exist in 
both attenuation and propagation zones of PE crystal. In spite approximate analytical solution for breathers 
has been obtained using a model of isolated chain, we confirmed numerically not only validity of such a 
model but revealed also that the parameters of breathers change unnoticeably if taking into account inter- 
chain interaction. The reason is a weakness of interchain interaction in comparison to intrachain one. Such 
a weakness is especially clear for optic excitations. The breathers demonstrated stability to mutual collision 
and have large enough time of life in the presence of thermal excitations. The interaction of optic breathers 
with supersonic and subsonic (topological) solitons may be both elastic and inelastic dependent on their 
parameters. 

The authors (L. I. M. and A. V. S.) thank the Region Rhone- Alpes and the Russian Foundation of Basic 
Research (awards 04-02-17306 and 04-03-321 19), respectively for financial support. 
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Figure 1: Schematic presentation of monoclinic structure of crystalline PE. The considered chain (curve 0) 
with local coordinates and 6 neighbor chains (curves 1-6) are shown. 
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Figure 2: Dispersion curves: optic (a, b) and acoustic (c, d) branches for S = 0.019 (a, c) and for S = 0.078 
(b, d) (parameter C 2 = 8C{). Thick part of the dispersion curves correspond to possible existence of the 
breathers. 
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Figure 3: Optical breather in PE chain (parameter 7 = 529 kJ/mol). The distribution along the chain energy 
E n (a), longitudinal w n+ i — w n (b) and transversal v n+ i — v n (c) relative displacements of the chain segments 
is shown (breather energy E = 56.06 kJ/mol, frequency u = 926 cm' 1 , velocity s = 0.15 (1146 m/s). 
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Figure 4: Optical breather in PE chain (parameter 7 = 130.122 kJ/mol). The distribution along the chain 
energy E n (a), longitudinal w n+ i — w n (b) and transversal v n+1 — v n (c) relative displacements of the chain 
segments is shown (breather energy E = 31.1 kJ/mol, frequency u = 909 cm -1 , velocity s = 0.48 (3888 
m/s). 
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Figure 5: Optical breather in PE chain (parameter 7 = 130.122 kJ/mol). The distribution along the chain 
energy E n (a), longitudinal w n+ \ — w n (b) and transversal v n+ i — v n (c) relative displacements of the chain 
segments is shown (breather energy E = 28.77 kJ/mol, frequency uj = 830 cm -1 , velocity s — 0.16 (1296 
m/s). 
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Figure 6: Periodic change of transversal u n and relative longitudinal displacements of zigzag w n+ \ — w n 
in the localization region of optic breathers under parameters (56), N = 200. Frequency of the breather 
uj = 820.5 cm" 1 is slightly lower than gap frequency. 
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Figure 7: Periodic change of transversal u n and relative longitudinal displacements of zigzag chain w n+ \ — 
w n in the localization region of optic breathers under parameter (57), N = 500. Frequency of breather 
uj = 928.4 cm -1 is slightly lower than frequency corresponding to low boundary of gap. 
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Figure 8: Breaking of the breather (frequency lu = 928.4 cm ) in thermalized chain (N = 500, T = 10 
K, A = 529 kJ/mol). Temporal dependence of current local magnitudes of temperature (kinetic energy of 
chain segments) T n is presented. 
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Figure 9: Decreasing the breather energy E in thermalized chain with A = 130.1 kJ/mol, to = 820.5 
cm" 1 (a) and A = 529 kJ/mol, uj = 928 cm" 1 (b) under temperature T = 1, 2, 3, 5, 10, 20 K (curves 1, 
2,. ..,6). Solid lines determine temporal dependencies of breather energy for concrete realization of chain 
thermalization, dotted lines - corresponding exponential law E(t) = E(0) exp(— at). 
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Figure 10: Collision of supersonic acoustic soliton with stationary optic breather (lu = 820.5 cm~l) (a), 
collision of moving breather (velocity s = 0.106s ) with stationary one (b) and elastic collision of two 
moving breathers (c). The temporal dependence of energy distribution E n along the chain is shown. 
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Figure 1 1 : Decreasing the breather energy E (a) and dimensionless breather velocity s/s (b) in PE chain. 
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Figure 12: Breaking of moving breather (velocity s = 0.106s ) in thermalized chain (temperature T = 5K). 
Temporal dependence of energy distribution along the chain E n is shown. 



29 




Figure 13: Passage of topological solitons (velocity s = 0.25s ) with charge q = (1, 0) (a), q = (0.5, 0.5) 
(b), and q = (0, 1) (c) through static breather (frequency u = 820.5 cm -1 ). 
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Figure 14: Breaking of the stationary breather (frequency uj = 820.5 cm' 1 ) as a result of collision with 
topological soli tons with integer charge q = (— 1, 0) (a) and half-integer charge q = (— 0.5, 0.5) (b) (veloc- 
ity of topological soliton s = 0.25s ). 
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